Language codes (from WALS) – WALS, ISO, ascii-name. These are used to merge across datasets.
codes = read.csv("../data/language_codes.csv") %>%
select(language, WALS, ISO) %>%
mutate(ISO = unlist(lapply(strsplit(as.character(ISO),
","),function(x) x[1])))
# in cases where there are two ISO codes, takes first
Lupyan & Dale (2010): Demographic variables and syntactic compelxity
wichmann_ISOS = read.csv("../data/ISO_mappings_from_wichmann.csv",
stringsAsFactors=F)[-1] %>%
unique()
ld = read.table("../data/lupyan_2010.txt", fill = T,
header = T, sep = "\t", na.strings = "*") %>%
left_join(codes, c("walsCode" = "WALS")) %>%
mutate(lang = tolower(lang)) %>%
left_join(wichmann_ISOS, by = "lang") %>% # use wichmann ISOS for some missing
mutate(ISO = as.factor(ifelse(is.na(ISO.x), ISO.y, ISO.x))) %>%
select(-ISO.y, -ISO.x) %>%
filter(!is.na(ISO)) %>%
filter(ISO != "")
# rename columns
ld = rename(ld, pop_log = logpop2,
mean.temp = aveTemp,
n.neighbors = numNeighbors,
sum.precip = sumPrecip,
sd.temp = sdTemp,
sd.precip = sdPrecip,
lang.family = langFamily,
lang.genus = langGenus,
native.country = nativeCountry,
native.country.area = nativeCountryArea,
growing.season = growingSeason)
# get means across language
ld.demo = ld %>%
group_by(ISO) %>%
summarise_each(funs(mean(., na.rm = TRUE)),
c(8:9, 16:121)) %>%
select(1:3, 93:95, 100:101, 103:105, 108)
# add in data family
fams = ld %>%
group_by(ISO) %>%
filter(row_number() == 1) %>%
select(ISO, lang.family, lang.genus, native.country,
native.country.area, language)
ld.demo = left_join(fams, ld.demo, by = "ISO") %>% ungroup()
WALS features of complexity
# get variables to include
qual = ld %>%
select(18:103, 124) %>%
mutate_each(funs(as.factor)) %>%
group_by(ISO) %>%
summarise_each(funs(most.frequent.level)) %>%
mutate_each(funs(as.factor))
ld_feature_names = read.csv("../data/lupyan_2010_all_feature_mappings.csv")
qualVarNames = intersect(names(qual), ld_feature_names$WALS.feature.name)
qual = select(qual,which(is.element(names(qual), c("ISO", qualVarNames))))
# remap factor levels to complexity values (0,1)
# note two variables that are reported in the paper (NICPOC[59] and CYSIND[39] are missing in the data)
for (i in 1:length(qualVarNames)){
thisVarLabs = ld_feature_names[ld_feature_names$WALS.feature.name == qualVarNames[i],]
old = thisVarLabs$ld.level.label
if (is.na(old)) {print("NA!!!")}
new = thisVarLabs$ld.complexity
col_i = grep(qualVarNames[i], colnames(qual))
qual[,col_i] = mapvalues(as.matrix(qual[,col_i]),
from = as.character(old),
to = as.character(new), warn_missing = TRUE)
}
ld.complexity = qual %>%
mutate_each(funs(as.numeric), -ISO) %>%
gather(variable, complexity.level, 2:28) %>%
group_by(ISO) %>%
summarise(morphological.complexity = sum(complexity.level))
ld.demo.qual = left_join(ld.demo, ld.complexity)
Bentz, et al. (2015): L2 learners and lexical diversity.Note here that we are only looking at languages from the UDHR corpus. The other corpora have smaller languages only.
bentz = read.csv("../data/bentz_2015.csv") %>%
gather(temp, LDT, starts_with("LDT")) %>%
unite(temp1, measure, temp, sep = ".") %>%
spread(temp1, LDT) %>%
filter(text == "UDHR") %>%
select(iso_639_3, RatioL2, a.LDTscaled,
H.LDTscaled, TTR.LDTscaled, Stock,
Region, L1_speakers, L2_speakers, TTR.LDT) %>%
rename(ISO = iso_639_3,
n.L1.speakers = L1_speakers,
n.L2.speakers = L2_speakers,
ratio.L2.L1 = RatioL2,
stock = Stock,
region = Region,
scaled.LDT.TTR = TTR.LDTscaled,
scaled.LDT.ZM = a.LDTscaled,
scaled.LDT.H = H.LDTscaled)
# n_speakers taken directly from ethnologue website by ML for missing languages that were high-frequency across the datasets
ethno_sup = read.csv("../data/supplementary_enthnologue.csv")
bentz = full_join(bentz, ethno_sup) %>%
mutate(ISO = as.factor(ISO))
Atkinson (2011)
atkinson = read.csv("../data/atkinson_2011.csv") %>%
select(normalized.vowel.diversity,
normalized.consonant.diversity,
normalized.tone.diversity,
normalized.phoneme.diversity, ISO,
distance.from.origin) %>%
mutate(ISO = as.factor(unlist(lapply(strsplit(as.character(ISO),
" "),function(x) x[1])))) %>%
distinct(ISO)
Lewis & Frank (under review): Complexity Bias
# complexity bias
cb = read.csv("../data/lewis_2015.csv") %>%
rename(complexity.bias = corr,
p.complexity.bias = p.corr,
mono.complexity.bias = mono.cor,
open.complexity.bias = open.cor) %>%
left_join(codes, by="language") %>%
select(-X.1, -X, -lower.ci, -upper.ci, -checked,
-mean.length) %>%
distinct(ISO) %>%
filter(language != "english") %>% # english is an outlier in this dataset because norms were colelcted in english
select(-language, -WALS)
Futrell, Mahowald, & Gibson (2015): Dependency length
dl = read.csv("../data/futrell_2015.csv") %>%
left_join(codes, by="language") %>%
select(-language, -WALS, -fixed.random.baseline.slope,
-observed.slope, -m_rand) %>%
rename(mean.dependency.length = m_obs) %>%
filter(!is.na(ISO))
Pellegrino, Coupe, & Marsico (2015): Information density
uid = read.csv("../data/pellegrino_2015.csv") %>%
left_join(codes, by="language") %>%
select(-language, -WALS)
Moran, McCloy, & Wright (2012): Phonemic inventory
phoneme = read.csv("../data/moran_2012.csv") %>%
select(ISO, pho, con, vow, son,
obs, mon, qua, ton) %>%
rename(n.phonemes = pho,
n.consonants = con,
n.vowels = vow,
n.sonorants = son,
n.obstruents = obs,
n.monophthongs = mon,
n.qual.monophthongs = qua,
n.tones = ton)
Wichmann, et al. (2013): Mean word length
# note is this IPA? - fix this!
# [308 low frequency language have no ISO, leaving 4421 (TOTAL:4729)]
ml.raw = read.csv("../data/wichmann_2013.csv") %>%
select(1,1:109) %>%
gather(word,translation,I:name) %>%
mutate(nchar = unlist(
lapply(
lapply(
strsplit(
gsub("[[:space:]]", "", translation) ,
","),
nchar), mean))) %>%
filter(translation != "")
# subset to only those words in the swadesh list (n = 40)
swadesh.words = ml[ml$ISO == "gwj", "word"]
ml.raw = ml.raw %>%
filter(is.element(word, swadesh.words))
ml = ml.raw %>%
group_by(ISO) %>%
summarize(mean.length = mean(nchar, na.rm = T)) %>%
filter(ISO != "")
#more_ISOs = ml.raw %>%
# select(ISO, names) %>%
# mutate(names = tolower(names))
Luniewska, et al. (2015): Age of acquistion (Aoa)
aoa.raw = read.csv("../data/luniewska_2015.csv",
header = T, fileEncoding = "latin1")
aoa = aoa.raw %>%
gather(language_aoa, aoa, grep("aoa",
names(aoa.raw))) %>%
mutate(language_aoa = tolower(unlist(lapply(strsplit(
as.character(language_aoa),"_"), function(x) x[1])))) %>%
select(-3:-27) %>%
rename(language = language_aoa) %>%
left_join(codes, by="language") %>%
mutate(language = as.factor(language)) %>%
group_by(ISO) %>%
summarize(mean.aoa = mean(aoa, na.rm = T)) %>%
filter(!is.na(ISO))
Merge data frames together by ISO
d = full_join(ld.demo.qual, cb, by = "ISO") %>%
full_join(bentz, by = "ISO") %>%
full_join(phoneme, by = "ISO") %>%
full_join(ml, by = "ISO") %>%
full_join(dl, by = "ISO") %>%
full_join(uid, by = "ISO") %>%
full_join(aoa, by = "ISO") %>%
full_join(atkinson, by = "ISO") %>%
mutate(ISO = as.factor(ISO)) %>%
filter(!is.na(native.country)) # this is important to have for the mixed effects models
load data
# write.csv(d, "../data/langLearnVar_data.csv")
d = read.csv("../data/langLearnVar_data.csv")[,-1]
Remove outliers and check for normality
d.clean = d %>%
mutate_each(funs(removeOutlier(.,N_EXCLUDE_SDS)), c(-ISO, -lang.family,
-native.country, -native.country.area,
-lang.genus, -language, -stock, -region))
d.clean %>%
select(-ISO, -lang.family, -native.country, -native.country.area,
-lang.genus, -language, -stock, -region) %>%
gather(var, value) %>%
ggplot(aes(sample = value)) +
stat_qq(na.rm = T, size = .2) +
facet_wrap( ~ var, scales = "free") +
theme_bw()
The following variables look right-skewed.
NN_vars = c("area", "perimeter", "n.neighbors", "sd.precip",
"ratio.L2.L1", "n.L1.speakers", "n.L2.speakers",
"n.phonemes", "n.consonants","n.vowels", "n.sonorants",
"n.obstruents", "n.monophthongs", "n.qual.monophthongs",
"n.tones")
Take the log of these variables, and check for normality again.
d.clean = d.clean %>%
mutate_each(funs(log = log(.), dummy = sum(.)),
one_of(NN_vars)) %>%
select(-contains("dummy")) %>%
select(-one_of(NN_vars))
d.clean %>%
select(-ISO, -lang.family, -native.country, -native.country.area,
-lang.genus, -language, -stock, -region) %>%
gather(var, value) %>%
ggplot(aes(sample = value)) +
stat_qq(na.rm = T, size = .2) +
facet_wrap( ~ var, scales = "free") +
theme_bw()
Check for normality of variables using Shapiro-Wilk Normality Test
is.na(d.clean) <- sapply(d.clean, is.infinite)
normality.test = d.clean %>%
summarise_each(funs(unlist(tidy(shapiro.test(.))[2])),
c(-ISO, -lang.family, -native.country,
-native.country.area,
-lang.genus, -language, -stock, -region)) %>%
gather(variable, shapiro.p.value) %>%
mutate(normal = ifelse(shapiro.p.value > .05, TRUE, FALSE)) %>%
summary()
By this test, only 11 variables normal. But, I think this is the best we can do (?).
Look at histograms of all variables.
demo_vars = c("lat", "lon", "perimeter_log", "n.neighbors_log",
"mean.temp", "sum.precip",
"sd.temp", "growing.season", "pop_log", "area_log",
"n.neighbors_log", "sd.precip_log", "ratio.L2.L1_log",
"n.L1.speakers_log",
"n.L2.speakers_log", "distance.from.origin")
lang_vars = c(setdiff(names(d.clean), demo_vars))[-c(1:6, 22,23)]
d.clean %>%
gather(variable, num, lat:length(d.clean)) %>%
filter(variable != "stock", variable != "region") %>%
mutate(var_type = ifelse(is.element(variable, demo_vars),
"demo", "lang"),
num = as.numeric(as.character(num))) %>%
arrange(var_type) %>%
mutate(variable = factor(variable, variable)) %>%
ggplot(aes(x = num, fill = var_type)) +
geom_histogram(position = "identity") +
facet_wrap( ~ variable, scales = "free") +
theme_bw()
Number of languages we have data, for each variable
d.clean %>%
gather(variable, num, lat:length(d.clean)) %>%
filter(variable != "stock", variable != "region") %>%
group_by(variable) %>%
summarize(counts = length(which(is.na(num) == F))) %>%
mutate(var_type = ifelse(is.element(variable, demo_vars),
"demo", "lang")) %>%
mutate(counts_trunc = ifelse(counts > 150, 150, counts)) %>%
arrange(var_type) %>%
mutate(variable = factor(variable, variable)) %>%
ggplot(aes(y=counts_trunc, x = variable, fill = var_type)) +
geom_bar(stat = "identity") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, hjust = 1),
legend.position = "none") +
ylab("Number of languages (truncated at 150)") +
ggtitle("Number languages per measure")
Look at number of languages intersecting each variable
arealControl_vars = c("stock", "region", "lang.family", "lang.genus",
"native.country", "native.country.area")
d.clean %>%
gather(lang_measure, lang_value,
which(is.element(names(.), c(lang_vars, arealControl_vars)))) %>%
group_by(lang_measure) %>%
gather(pop_measure, pop_value,
which(is.element(names(.), demo_vars))) %>%
group_by(lang_measure, pop_measure) %>%
mutate(both_not_na = !is.na(lang_value) & !is.na(pop_value)) %>%
summarise(n_languages = length(which(both_not_na == TRUE))) %>%
ggplot(aes(pop_measure, lang_measure, group = lang_measure)) +
geom_tile(aes(fill = n_languages)) +
geom_text(aes(fill = n_languages, label = n_languages)) +
scale_fill_gradient(low = "white", high = "red") +
xlab("population variables") +
ylab("language variables") +
ggtitle("Number languages between variables") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
Fits controling for family in a mixed effect model
# write.csv(d.clean,"../data/langLearnVar_data_clean.csv")
# d.clean = read.csv("../data/langLearnVar_data_clean.csv")[-1]
# these are the variables we care about for the paper
demo_vars_crit = c("n.neighbors_log",
"mean.temp", "sum.precip", "sd.temp",
"pop_log", "area_log", "n.neighbors_log",
"sd.precip_log", "ratio.L2.L1_log")
lang_vars_crit = c("p.complexity.bias","TTR.LDT",
"mean.dependency.length", "mean.length",
"n.consonants_log", "n.vowels_log",
"morphological.complexity")
d.scatter.fam = d.clean %>%
select(pop_log, ratio.L2.L1_log,
n.neighbors_log, area_log, mean.temp, sd.temp, sum.precip,
sd.precip_log, morphological.complexity,TTR.LDT, p.complexity.bias,
mean.length,n.consonants_log,n.vowels_log,
lang.family, native.country) %>%
#mutate_each(funs(as.numeric(scale(.))), c(-15,-16)) %>% #scaling doesn't change anything
gather(lang_measure, lang_value,
which(is.element(names(.), lang_vars_crit))) %>%
group_by(lang_measure) %>%
gather(pop_measure, pop_value,
which(is.element(names(.), demo_vars_crit)))
is.na(d.scatter.fam) <- sapply(d.scatter.fam, is.infinite)
d.scatter.fam <- filter(d.scatter.fam, !is.na(pop_value),
!is.na(lang_value))
# get model fits
d.model.fits = d.scatter.fam %>%
group_by(lang_measure, pop_measure) %>%
do(tidy(lmer(lang_value ~ pop_value + (pop_value|lang.family) +
(1|native.country), data=.)))
line.fits = d.model.fits %>%
filter(term == "pop_value" | term == "(Intercept)") %>%
select(-std.error, -statistic, -group) %>%
spread(term, estimate) %>%
rename(s = pop_value)
names(line.fits)[3] = "i"
sigs = d.model.fits %>%
filter(term == "pop_value") %>%
mutate(sig.col = ifelse(statistic > 1.96, "pos",
ifelse(statistic < -1.96, "neg",
"none"))) %>%
mutate(pop_value = .1, lang_value = .1) %>% # this is a hack
ungroup
Paper plot #1
facet_pop_names <- c(
pop_log="Population \n(log)",
ratio.L2.L1_log="L2:L1 \n(log)",
n.neighbors_log="N. Ling. \nNeighbors (log)",
area_log="Area \n(log)",
mean.temp="Temperature \n Mean",
sd.temp="Temperature \nSD",
sum.precip="Precipitation \nSum",
sd.precip_log="Precipitation \nSD (log)"
)
facet_lang_names <- c(
morphological.complexity="Morphosyntactic \nComplexity",
TTR.LDT="Lexical \nDiversity",
p.complexity.bias="Complexity \nBias",
mean.length="Word \nLength",
n.consonants_log="N. Consonants \n(log)",
n.vowels_log="N. Vowels \n(log)"
)
pdf("../papers/cogsci2016/figs/plot1.pdf", width = 18, height = 14)
ggplot(d.scatter.fam, aes(x = pop_value, y = lang_value)) +
geom_rect(data = sigs, aes(fill = sig.col),
xmin = -Inf, xmax = Inf,
ymin = -Inf, ymax = Inf, alpha = 0.2) +
geom_point(size = .3) +
geom_abline(line.fits, mapping = aes(slope = s, intercept = i),
size = 2, color = "green") +
facet_grid(lang_measure ~ pop_measure, scales = "free",
labeller = labeller(pop_measure = facet_pop_names,
lang_measure = facet_lang_names)) +
scale_fill_manual(values = c( "mediumblue", "grey99", "red1")) +
xlab("\nEnvironmental variables") +
ylab("Language variables\n") +
theme_bw() +
theme(legend.position = "none",
axis.title = element_text(size = 21),
strip.text= element_text(size = 17))
dev.off()
prep by checking for correlations
# get demographic variables only (excluding L2 because intersections contains only 40 languages)
d.demo = d.clean %>%
select(lat, mean.temp, sd.temp, sum.precip, sd.precip_log,
area_log, pop_log, n.neighbors_log) %>%
mutate(lat = abs(lat))
is.na(d.demo) <- sapply(d.demo, is.infinite)
demo.corr = cor(d.demo, use = "pairwise.complete")
abs(demo.corr) > CORR_THRESHOLD
## lat mean.temp sd.temp sum.precip sd.precip_log area_log
## lat TRUE TRUE FALSE FALSE TRUE FALSE
## mean.temp TRUE TRUE FALSE FALSE FALSE FALSE
## sd.temp FALSE FALSE TRUE FALSE FALSE FALSE
## sum.precip FALSE FALSE FALSE TRUE FALSE FALSE
## sd.precip_log TRUE FALSE FALSE FALSE TRUE FALSE
## area_log FALSE FALSE FALSE FALSE FALSE TRUE
## pop_log FALSE FALSE FALSE FALSE FALSE FALSE
## n.neighbors_log FALSE FALSE FALSE FALSE FALSE FALSE
## pop_log n.neighbors_log
## lat FALSE FALSE
## mean.temp FALSE FALSE
## sd.temp FALSE FALSE
## sum.precip FALSE FALSE
## sd.precip_log FALSE FALSE
## area_log FALSE FALSE
## pop_log TRUE FALSE
## n.neighbors_log FALSE TRUE
# exclude lat (correlated with mean.temp, sd.precip_log) [threshold: > .8]
PCA
d.demo = select(d.demo, -lat)
# do pca
pca.demo = prcomp(na.omit(d.demo), scale = TRUE)
# look at variance explained
fviz_eig(pca.demo, addlabels = TRUE,
linecolor ="red", geom = "line") +
theme_bw()
barplot(pca.demo$sdev/pca.demo$sdev[1])
pca.demo = prcomp(na.omit(d.demo), scale = TRUE, tol = .6)
# plot loadings
pcs.demo <- cbind(names(d.demo),
gather(as.data.frame(pca.demo$rotation),
pc, value))
names(pcs.demo)[1] = "variable"
pcs.demo$variable = factor(pcs.demo$variable, levels = pcs.demo$variable[1:8])
# order variables
pcs.demo$magnitude = ifelse(abs(pcs.demo$value)>.2, TRUE, FALSE)
ggplot(pcs.demo) +
geom_bar(aes(x = variable, y = value,
fill = magnitude), stat="identity") +
facet_wrap(~pc) +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
scale_fill_manual(name = "Magnitude",
values = c("grey","red1")) +
ggtitle("PCA loadings")
# merge PCAs with full dataset for model fits
d.demo[!is.na(d.demo)] = 0
d.demo[is.na(d.demo)] = 1
good_is = which(rowSums(d.demo) == 0)
d.pca.demo = d.clean[good_is,] %>% cbind(pca.demo$x)
# flip size PC so can be interpretted as "bigness"
#d.pca.demo$PC2 = -d.pca.demo$PC2
# rename PCs
names(d.pca.demo)[names(d.pca.demo) == "PC1"] = "hot"
names(d.pca.demo)[names(d.pca.demo) == "PC2"] = "small"
Maps of principle components for demographic variable by language
ggplot(d.pca.demo) +
borders("world", colour="gray50", fill="gray50") +
geom_point(aes(x = lon, y=lat,
color = hot), size=3) +
scale_colour_gradient(high="green", low = "white") +
ggtitle("Hot (PC1)") +
mapTheme
ggplot(d.pca.demo) +
borders("world", colour="gray50", fill="gray50") +
geom_point(aes(x = lon, y=lat,
color = small), size=3) +
scale_colour_gradient(high="blue", low = "white") +
ggtitle("Small (PC2)") +
mapTheme
Model fits using demographic PCA with all language variables
demo_vars_crit2 = c("hot", "small")
lang_vars_crit2 = c("morphological.complexity",
"p.complexity.bias","TTR.LDT",
"mean.length", "mean.aoa",
"n.consonants_log", "n.vowels_log")
d.scatter.fam = d.pca.demo %>%
select(n.consonants_log, n.vowels_log,
mean.length, p.complexity.bias, TTR.LDT,
morphological.complexity,
lang.family, native.country, small, hot) %>%
gather(lang_measure, lang_value,
which(is.element(names(.), lang_vars_crit2))) %>%
group_by(lang_measure) %>%
gather(pop_measure, pop_value,
which(is.element(names(.), demo_vars_crit2)))
is.na(d.scatter.fam) <- sapply(d.scatter.fam, is.infinite)
d.scatter.fam <- filter(d.scatter.fam, !is.na(pop_value),
!is.na(lang_value))
# get model fits
d.model.fits = d.scatter.fam %>%
group_by(lang_measure, pop_measure) %>%
do(tidy(lmer(lang_value ~ pop_value +
(pop_value|lang.family) +
(1|native.country),
data=.)))
line.fits = d.model.fits %>%
filter(term == "pop_value" | term == "(Intercept)") %>%
select(-std.error, -statistic, -group) %>%
spread(term, estimate) %>%
rename(s = pop_value)
names(line.fits)[3] = "i"
sigs = d.model.fits %>%
filter(term == "pop_value") %>%
mutate(sig.col = ifelse(statistic > 1.96, "pos",
ifelse(statistic < -1.96, "neg",
"none"))) %>%
mutate(pop_value = .1, lang_value = .1) %>% # this is a hack
ungroup
ggplot(d.scatter.fam, aes(x = pop_value, y = lang_value)) +
geom_rect(data = sigs, aes(fill = sig.col),
xmin = -Inf, xmax = Inf,
ymin = -Inf, ymax = Inf, alpha = 0.2) +
geom_point(size = .3) +
geom_abline(line.fits, mapping = aes(slope = s, intercept = i), size = 1, color = "green") +
facet_grid(lang_measure~pop_measure, scales = "free") +
scale_fill_manual(values = c( "mediumblue", "grey99","red1")) +
theme_bw() +
xlab("Demographic variables") +
ylab("Language variables") +
theme(legend.position = "none")
prep by looking at correlations
d.lang = d.clean %>%
select(morphological.complexity, scaled.LDT.H,
mean.length, n.consonants_log, n.vowels_log)
is.na(d.lang) <- sapply(d.lang, is.infinite)
lang.corr = cor(d.lang, use = "pairwise.complete")
abs(lang.corr) > CORR_THRESHOLD
## morphological.complexity scaled.LDT.H mean.length
## morphological.complexity TRUE FALSE FALSE
## scaled.LDT.H FALSE TRUE FALSE
## mean.length FALSE FALSE TRUE
## n.consonants_log FALSE FALSE FALSE
## n.vowels_log FALSE FALSE FALSE
## n.consonants_log n.vowels_log
## morphological.complexity FALSE FALSE
## scaled.LDT.H FALSE FALSE
## mean.length FALSE FALSE
## n.consonants_log TRUE FALSE
## n.vowels_log FALSE TRUE
# include all language variables
PCA
# do pca
pca.lang = prcomp(na.omit(d.lang), scale = TRUE)
barplot(pca.lang$sdev/pca.lang$sdev[1])
# look at variance explained
fviz_eig(pca.lang, addlabels = TRUE,
linecolor ="red", geom = "line") +
theme_bw()
pca.lang = prcomp(na.omit(d.lang), scale = TRUE, tol = .7)
# plot loadings
pcs.lang <- cbind(names(d.lang),
gather(as.data.frame(pca.lang$rotation),
pc, value))
names(pcs.lang)[1] = "variable"
pcs.lang$variable = factor(pcs.lang$variable , levels = pcs.lang$variable[1:8]) # order variables
pcs.lang$magnitude = ifelse(abs(pcs.lang$value)>.2, TRUE, FALSE)
ggplot(pcs.lang) +
geom_bar(aes(x = variable, y = value,
fill = magnitude), stat = "identity") +
facet_wrap(~pc) +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
scale_fill_manual(name = "Magnitude",
values = c("grey","red1")) +
ggtitle("PCA loadings")
# merge PCAs with full dataset for model fits
d.lang[!is.na(d.lang)] = 0
d.lang[is.na(d.lang)] = 1
good_is = which(rowSums(d.lang) == 0)
d.pca.lang = d.clean[good_is,] %>% cbind(pca.lang$x)
names(d.pca.lang)[names(d.pca.lang) == "PC1"] = "PC1_L"
names(d.pca.lang)[names(d.pca.lang) == "PC2"] = "PC2_L"
Model fits using language PCA with all demographics variables
lang_vars_crit2 = c("PC1_L", "PC2_L")
demo_vars_crit2 = c("n.neighbors_log",
"mean.temp", "sum.precip", "sd.temp",
"pop_log", "area_log", "n.neighbors_log",
"sd.precip_log", "ratio.L2.L1_log")
d.scatter.fam = d.pca.lang %>%
select(ratio.L2.L1_log, pop_log,
n.neighbors_log, area_log, mean.temp, sd.temp, sum.precip,
PC1_L, PC2_L, lang.family, native.country) %>%
gather(lang_measure, lang_value,
which(is.element(names(.), lang_vars_crit2))) %>%
group_by(lang_measure) %>%
gather(pop_measure, pop_value,
which(is.element(names(.), demo_vars_crit2)))
is.na(d.scatter.fam) <- sapply(d.scatter.fam, is.infinite)
d.scatter.fam <- filter(d.scatter.fam, !is.na(pop_value),
!is.na(lang_value))
# get model fits
d.model.fits = d.scatter.fam %>%
group_by(lang_measure, pop_measure) %>%
do(tidy(lmer(lang_value ~ pop_value +
(1|native.country), data=.)))
line.fits = d.model.fits %>%
filter(term == "pop_value" | term == "(Intercept)") %>%
select(-std.error, -statistic, -group) %>%
spread(term, estimate) %>%
rename(s = pop_value)
names(line.fits)[3] = "i"
sigs = d.model.fits %>%
filter(term == "pop_value") %>%
mutate(sig.col = ifelse(statistic > 1.96, "pos",
ifelse(statistic < -1.96, "neg",
"none"))) %>%
mutate(pop_value = .1, lang_value = .1) %>% # this is a hack
ungroup
ggplot(d.scatter.fam, aes(x = pop_value, y = lang_value)) +
geom_rect(data = sigs, aes(fill = sig.col),
xmin = -Inf, xmax = Inf,
ymin = -Inf, ymax = Inf, alpha = 0.2) +
geom_point(size = .3) +
geom_abline(line.fits, mapping = aes(slope = s, intercept = i), size = 1, color = "green") +
facet_grid(lang_measure~pop_measure, scales = "free") +
scale_fill_manual(values = c( "mediumblue", "grey99","red1")) +
theme_bw() +
xlab("Demographic variables") +
ylab("Language variables") +
theme(legend.position = "none")
Model fits
all.pca = left_join(d.pca.demo,
d.pca.lang[,c("PC1_L", "PC2_L", "ISO")],
by = "ISO")
demo_vars_crit2 = c("hot", "small")
lang_vars_crit2 = c("PC1_L", "PC2_L")
d.scatter.fam = all.pca %>%
select(PC1_L, PC2_L, hot, small,
lang.family, native.country, ISO, language) %>%
gather(lang_measure, lang_value,
which(is.element(names(.), lang_vars_crit2))) %>%
group_by(lang_measure) %>%
gather(pop_measure, pop_value,
which(is.element(names(.), demo_vars_crit2)))
is.na(d.scatter.fam) <- sapply(d.scatter.fam, is.infinite)
d.scatter.fam <- filter(d.scatter.fam, !is.na(pop_value),
!is.na(lang_value))
# get model fits
d.model.fits = d.scatter.fam %>%
group_by(lang_measure, pop_measure) %>%
do(tidy(lmer(lang_value ~ pop_value + (pop_value|lang.family) +
(1|native.country), data=.)))
line.fits = d.model.fits %>%
filter(term == "pop_value" | term == "(Intercept)") %>%
select(-std.error, -statistic, -group) %>%
spread(term, estimate) %>%
rename(s = pop_value)
names(line.fits)[3] = "i"
sigs = d.model.fits %>%
filter(term == "pop_value") %>%
mutate(sig.col = ifelse(statistic > 1.96, "pos",
ifelse(statistic < -1.96, "neg",
"none"))) %>%
mutate(pop_value = .1, lang_value = .1) %>% # this is a hack
ungroup
ggplot(d.scatter.fam, aes(x = pop_value, y = lang_value)) +
geom_rect(data = sigs, aes(fill = sig.col),
xmin = -Inf, xmax = Inf,
ymin = -Inf, ymax = Inf, alpha = 0.2) +
geom_point(size = .3) +
geom_abline(line.fits, mapping = aes(slope = s, intercept = i), size = 1, color = "green") +
facet_grid(lang_measure~pop_measure, scales = "free") +
scale_fill_manual(values = c( "mediumblue", "grey99","red1")) +
theme_bw() +
xlab("Demographic variables") +
ylab("Language variables") +
theme(legend.position = "none")
Paper plot #2
vplayout <- function(x, y)
viewport(layout.pos.row = x, layout.pos.col = y)
facet_env_names <- c(PC1="PC1: Hot and rainy", PC2="PC2: Small")
facet_lang_names <- c(PC1="PC1: Complex", PC2="PC2: Short words")
pc.demo = ggplot(pcs.demo) +
geom_bar(aes(x = variable, y = value, fill = variable), stat="identity") +
facet_wrap(~pc, labeller = labeller(pc = facet_env_names)) +
ylab("Component loading")+
theme_bw() +
theme(axis.text.x = element_blank(),
axis.title.x = element_blank(),
strip.text= element_text(size = 21),
axis.title = element_text(size = 21)) +
scale_x_discrete(breaks=NULL) +
scale_fill_brewer(palette="YlOrRd") +
#theme(legend.margin=unit(-5,"cm"))
theme(legend.position = "none")
pc.lang = ggplot(pcs.lang) +
geom_bar(aes(x = variable, y = value,
fill = variable), stat="identity") +
theme_bw() +
theme(axis.text.y = element_blank(),
axis.title.y = element_blank(),
strip.text= element_text(size = 21),
axis.title = element_text(size = 21)) +
coord_flip() +
facet_grid(pc ~ ., scales="free_x",
labeller = labeller(pc = facet_lang_names)) +
scale_x_discrete(breaks=NULL) +
ylab("Component loading") +
theme(legend.position = "none") +
scale_fill_brewer(palette="BuGn")
p1 <- ggplot(d.scatter.fam, aes(x = pop_value, y = lang_value)) +
geom_rect(data = sigs, aes(fill = sig.col),
xmin = -Inf, xmax = Inf,
ymin = -Inf, ymax = Inf, alpha = 0.2) +
geom_point() +
geom_abline(line.fits, mapping = aes(slope = s, intercept = i),
size = 2) +
facet_grid(lang_measure~pop_measure, scales = "free") +
geom_text_repel(aes(label = language)) +
scale_fill_manual(values = c( "mediumblue", "grey99","red1")) +
theme_bw() +
theme(strip.background = element_blank(),
strip.text = element_blank(),
axis.title = element_text(size = 21)) +
ylab("Rotated language values")+
xlab("Rotated environmental values")+
theme(legend.position = "none")
#pdf("../papers/cogsci2016/figs/plot2.pdf", width = 18, height = 12)
grid.newpage()
pushViewport(viewport(layout = grid.layout(6, 6)))
print(pc.demo, vp=vplayout(1:2.3,1:4))
print(pc.lang, vp=vplayout(3:6,5:6))
print(p1, vp=vplayout(3:6,1:4))
#dev.off()
# predicted relationship with AOA with TTR and number of consonants
# morphological complexity
tidy(lmer(mean.aoa ~ morphological.complexity + (morphological.complexity|lang.family) + (1|native.country), d.clean))[2, "statistic"]
cor.test(d.clean$mean.aoa, d.clean$morphological.complexity)
# complexity bias
tidy(lmer(mean.aoa ~ p.complexity.bias + (p.complexity.bias|lang.family) + (1|native.country), d.clean))[2, "statistic"]
cor.test(d.clean$mean.aoa, d.clean$p.complexity.bias)
# TTR
tidy(lmer(mean.aoa ~ TTR.LDT+ (TTR.LDT|lang.family) + (1|native.country), d.clean))[2, "statistic"]
cor.test(d.clean$mean.aoa, d.clean$scaled.LDT.TTR)
cor.test(d.clean$mean.aoa, d.clean$scaled.LDT.H)
cor.test(d.clean$mean.aoa, d.clean$scaled.LDT.ZM)
# dependency length
tidy(lmer(mean.aoa ~ mean.dependency.length + (mean.dependency.length|lang.family) + (1|native.country), d.clean))[2, "statistic"]
cor.test(d.clean$mean.aoa, d.clean$mean.dependency.length)
# mean length
tidy(lmer(mean.aoa ~ mean.length + (mean.length|lang.family) + (1|native.country), d.clean)) [2, "statistic"]
cor.test(d.clean$mean.aoa, d.clean$mean.length)
# n_consonants **
tidy(lmer(mean.aoa ~ n.consonants_log + (n.consonants_log|lang.family) + (1|native.country) , d.clean))[2, "statistic"]
tidy(lmer(mean.aoa ~ normalized.consonant.diversity + (1|native.country) , d.clean)) [2, "statistic"]
# n_phonemes
tidy(lmer(mean.aoa ~ n.phonemes_log + (n.phonemes_log|lang.family) + (1|native.country) , d.clean))[2, "statistic"]
#n_vowels
tidy(lmer(mean.aoa ~ n.vowels_log + (n.vowels_log|lang.family) + (1|native.country) , d.clean))[2, "statistic"]
ggplot(d.clean, aes(mean.length, mean.aoa, label = ISO)) +
geom_point() +
geom_label() +
geom_smooth(method = "lm")
ggplot(d.clean, aes(scaled.LDT.TTR, mean.aoa, label = ISO)) +
geom_point() +
geom_label() +
geom_smooth(method = "lm")
ggplot(d.clean, aes(n.consonants_log, mean.aoa, label = ISO)) +
geom_point() +
geom_label() +
geom_smooth(method = "lm")
ggplot(d.clean, aes(n.phonemes_log, mean.aoa, label = ISO)) +
geom_point() +
geom_label() +
geom_smooth(method = "lm")
ggplot(d.clean, aes(n.vowels_log, mean.aoa, label = ISO)) +
geom_point() +
geom_label() +
geom_smooth(method = "lm")
ggplot(d.clean, aes(mean.dependency.length , mean.aoa, label = ISO)) +
geom_point() +
geom_label() +
geom_smooth(method = "lm")
d.clean = ld %>%
select(ISO, DRYSOV) %>%
left_join(d.clean, by="ISO")
old = 1:7
new = new = c("SOV", "SVO", "VSO", "VOS", "OVS", "OSV", "none")
d.clean$DRYSOV = mapvalues(d.clean$DRYSOV,
from = as.character(old),
to = as.character(new))
d.clean %>%
distinct(ISO) %>%
filter(!is.na(mean.aoa))%>%
group_by(DRYSOV) %>%
multi_boot_standard("mean.aoa", na.rm = T) %>%
ggplot(aes(x = DRYSOV, y = mean, fill = DRYSOV)) +
geom_bar(stat = "identity") +
geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper),
size=.3,
width=.2,
position=position_dodge(.9))
ld.demo.qual %>%
mutate(morphological.complexity.factor =
as.factor(morphological.complexity)) %>%
group_by(morphological.complexity.factor) %>%
multi_boot(column="pop_log",
summary_groups = c("morphological.complexity.factor"),
statistics_functions = c("mean", "ci_lower","ci_upper")) %>%
mutate(morphological.complexity =
as.numeric(as.character(morphological.complexity.factor))) %>%
ggplot() +
geom_pointrange(aes(x= morphological.complexity.factor,
y = mean, ymin = ci_lower, ymax = ci_upper)) +
geom_smooth(method = "lm",aes(morphological.complexity-5, mean)) +
# theres a weird offset because complexity is a factor?
ylab("Log population") +
xlab("Morphological complexity") +
ggtitle("Lupyan and Dale (2010) Fig. 3 reproduction") +
theme_bw()